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ABSTRACT 


This  technical  memorandum  describes  a  stochastic  modeling 
process  that  will  enable  accurate  determination  of  the  misclassi- 
fication  probabilities  and  performance  loss  caused  by  use  of  a 
digital  multibeam  steering  system  (DIMUS)  instead  of  a  linear 
signal  processor  during  beamforming.  A  compact  program  for  the 
evaluation  of  the  basic  clipped  correlation  function  is  also 
derived  and  presented. 
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COMPARISON  OF  LINEAR  AND  DIMUS  BEAMFORMING 
IN  TERMS  OF  MISCLASSIFI CATION  PROBABILITIES 
FOR  A  SET  OF  DETERMINISTIC  SIGNALS 

INTRODUCTION 

An  idealized  communications  scenario  will  be  posed  here  in 
order  to  accurately  determine  the  performance  losses  caused  by 
the  use  of  a  digital  multibeam  steering  system  (DIMUS)  instead  of 
a  linear  signal  processor  during  beamforming.  These  relative 
losses  are  of  great  interest  in  communications,  especially  for 
moderate-to-high  input  signal-to-noi se  ratios  (SNRs),  where 
saturation  and  distortion  effects  can  take  place  in  the  clippers. 
Although  the  absolute  levels  of  performance  achieved  by  the 
following  idealized  models  of  linear  and  DIMUS  beamformers 
(combined  with  crosscorrelators)  may  not  be  attainable  in 
practice,  the  difference  in  levels  should  be  an  accurate  measure 
of  the  loss  incurred  by  the  use  of  clipping.  This  difference  can 
be  investigated  as  a  function  of  the  input  SNR,  the  number  of 
channels,  the  number  of  time  samples,  the  reference  waveforms, 
etc.,  in  order  to  pinpoint  and  ascertain  where  further  effort  or 
modified  processing  may  be  required  to  minimize  the  effects  of 
the  additional  losses. 

A  set  of  M  deterministic  real  signals  (sm(t)}  is  used  to 
communicate  information  to  a  receiver,  which  knows  these  detailed 
waveforms,  but  not  which  one  was  transmitted.  If  signal  m  is 
transmitted,  the  received  waveform  at  the  k-th  element  of  the 
receiver  (after  beamformer  delay-time  alignment)  is 
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t 


(1) 


s  (t)  +  n,  (t)  for  1  <  k  <  K 
m  k  —  — 

where  additive  noises  {n^(t)},  1  <  k  <  K,  are  stationary,  real, 
zero-mean,  dependent  Gaussian  random  processes  with  covariances 

E{nk(t)  n j ( t-T) }  =  Rkj(T)  for  1  <  k,j  <  K  .  (2) 

The  quantity  E{  }  denotes  an  ensemble  average.  Based  upon  the  K 
measurements  of  the  waveforms  in  equation  (1)  over  an  observation 
time  T,  a  decision  must  be  made  as  to  which  of  the  M  signals 
{ sm ( t ) }  was  actually  transmitted. 

The  signal  set  { sm ( t ) }  could  be  augmented  with  an  additional 
signal,  sQ(t)  =  0,  in  order  to  include  the  detection  scenario, 
where  it  is  not  known  if  any  signal  is  present  at  all.  That 
situation  can  easily  be  incorporated  in  the  following 
developments,  if  desired,  by  a  slight  modification  to  the 
receiver  signal  processing. 
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LINEAR  PROCESSING 


BACKGROUND 

The  processor  of  interest  in  this  section  is  depicted  in 
figure  1.  Time  is  normalized  and  discretized  such  that  samples 
are  taken  at  multiples  of  unit  time,  at  instants  t  =  1,2,...,T. 


^(t) 


K  Inputs  Sum  on 
Space 


Multiply  Sum  on  M  Output 
References  Time  Variables 


Figure  1.  Linear  Processor 


The  m— th  local  reference,  r  (t),  would  normally  be  taken  as 
some  scaled  version  of  signal  s  (t),  at  least  if  the  received 
additive  noises  are  independent  with  white  spectra.  Generality 
will  be  allowed  here  by  not  specifying  the  exact  choices  of 
references  This  permits  reference  rm(t)  to  be  partially 

mismatched  to  the  m-th  signal,  such  as  a  linearly  filtered 
version  of  signal  sm(t);  a  clipped  version  of  the  signal,  such  as 
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sgn{sin(t)};  or  some  other  distorted  version,  such  as  g  { s  ( t ) } ,  if 
desired . 

The  output  variables  Z  s  [z^,...,zM]  in  figure  1  are  subject 
to  additional  processing  before  a  decision  is  made  as  to  which 
particular  signal  was  transmitted.  In  particular,  the  a  priori 
probabilities  of  each  signal,  the  energies  of  each  signal,  and 
the  costs  of  misclassif ication  of  each  signal  pair,  1  <  m,m  <  M, 
must  be  taken  into  account.  The  end  result  is  that  Mxl  output 
vector  Z  is  transformed  to  the  Mxl  decision  vector  Z  according  to 
linear  transformation 

Z  =  A  Z  +  B  ,  (3) 

where  known  MxM  matrix  A  and  known  Mxl  matrix  B  incorporate  the 
probabilities,  energies,  and  costs  mentioned  above.  The  signal 
identity  decision  then  consists  of  selection  of  the  number  of  the 
maximum  variable  in  modified  set  Z. 

The  noise  inputs  {nk(t)}  in  figure  1  are  joint  Gaussian 
random  processes.  Because  all  the  operations  on  the  inputs  in 
figure  1  are  linear,  the  output  variables  Z  are  also  joint 
Gaussian.  Therefore,  the  modified  decision  variables  Z  generated 
by  equation  (3)  are  joint  Gaussian,  with  mean  and  covariance 
matrices  given  by 

E{Z}  =  A  E { Z }  +  B  ,  Cov { Z }  =  A  Cov{Z}  AT  ,  (4) 
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respectively.  This  knowledge  in  equation  (4)  constitutes  a 
complete  statistical  description  of  the  decision  variables  Z. 

The  essential  analysis  problem  relevant  to  the  linear 
processor  in  figure  1  is  therefore  the  determination  of  the  mean 
and  covariance  matrices,  namely,  E{Z}  and  Cov{Z},  of  output  Z. 

SIMULATION 

Once  the  information  in  equation  (4)  is  available,  the 
probability  of  correct  decision  for  the  m-th  transmitted  signal, 
Pc(m),  requires  integration  of  the  joint  Gaussian  probability 
density  function  (PDF)  of  Z  over  a  sector  of  M-dimensional  space. 
This  is  an  intractable  analytic  problem,  except  in  the  most 
special  of  cases.  Accordingly,  it  is  recommended  that  the 
M-dimensional  vector  Z,  with  the  known  statistics  (4),  be 
simulated  for  purposes  of  estimating  P  (m) .  Of  course,  this 

V 

simulation  must  be  conducted  separately  for  each  value  of  m  in 
the  range  (1,M),  since  the  energies,  costs,  and  a  priori 
probabilities  will  be  different.  However,  this  is  not  a  major 
numerical  task  because  the  size  M  of  the  signal  set  is  less  than 
10,  and  the  correct-decision  probabilities  of  interest  are  on  the 
order  of  0.9  to  0.99.  This  contrasts  with  the  estimation  of 
false  alarm  probabilities,  which  can  be  in  the  10-6  range.  The 
generation  of  random  vector  Z  with  specified  statistics  (4)  is 
addressed  in  appendix  A. 
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This  limited  simulation  requirement  is  achieved  because  the 
linear  processor  in  figure  1  has  been  analyzed  exactly  up  through 
the  output  variables  Z.  Thus,  it  is  not  necessary  to  simulate  T 
time  samples  of  K  channels  of  correlated  input  data  waveforms, 
which  significantly  reduces  the  size  of  the  simulation.  Effort 
can  be  concentrated  instead  on  a  small-size  M-ary  decision,  where 
the  fraction  of  correct  decisions  for  signal  m  can  be  accurately 
estimated  without  an  excessive  simulation  schedule.  Integer  m 
must  be  varied  over  the  range  1  <  m  <  M. 


ANALYSIS 

The  output  variables  {zm}  in  figure  1  are  given  by  linear 
operations  on  the  inputs  according  to 

zm  -  E  zB<t>  TZ  [»«<*>  +  «k(t>]  - 

T  T 

=  K  YU  r  (t)  s  ( t )  +  YU  r_(t)  n<t)  for  1  <  rn  <  M  ,  (5) 

t=l  -  t=l 


where  n(t)  is  the  sum  of  all  K  received  noises  {n^(t)}.  The  mean 
of  output  zm  is  given  by  the  first  term  as 


T 

E { z  }  =  K  YZ  r  (t)  s( t ) 
m  t— 1  111  — 


for  1  <  m  <  M  , 


(6) 


namely,  K  times  the  crosscorrelation  of  reference  waveform  m  with 
transmitted  signal  m. 
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The  covariance  of  output  random  variables  z  and  z»  follows 

m  a 

from  the  second  term  in  equation  (5)  as 


E { zm  =  C  rm(t)  rX(u)  E{n(t>  n(u)J  = 

t  ,u=l 


=  JZ2  R(T)  for  1  ^  ^  M  • 


t=-T 


(7) 


where 


K 


K 


R(t)  2  E{n(t)  n(t-x)}  =  E{n.(t)  n.(t-x)}  =  YU  ^(t)  (8) 


k  ,  j=l 


r-1  .  k  j 

k ,3=1  J 


is  the  autocovariance  of  the  total  received  stationary  noise 
process  n(t),  and  where  equations  (5)  and  (2)  have  been  used 
The  quantity  is  the  autocorrelation  of  the  references 
according  to 


«J>mX('c)  "  H  rm{t)  r*(t“T)  for  1  ^  ^  M  ' 


m 


(9) 


the  sum  on  t  is  taken  wherever  the  summand  is  nonzero. 


Equations  (6)  and  (7)  furnish  the  mean  and  covariance 
information  required  to  statistically  characterize  output  random 
variables  { zm }  in  figure  1.  They  enable  the  simulation  described 
above  to  be  carried  out  for  any  signal  and  reference  sets.  All 
the  noise  covariances  in  equation  (2)  must  also  be  available; 
recall  that  these  noises  { ( t ) }  in  equation  (1)  and  figure  1  are 
measured  at  the  outputs  of  the  steering  delays  in  the  beamformer. 
In  terms  of  the  measured  element  noises  {n^(t)},  there  follows 
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nk(t)  =  ( t-Tjc(  9 )  ) ,  where  1^(0)  is  the  k-th  time  delay  required 

to  steer  in  direction  6. 


DIMUS  PROCESSING 


BACKGROUND 

In  this  section,  the  linear  processor  of  figure  1  is  modified 
to  include  clipping  at  the  element  level.  The  resulting  DIMUS 
processor  of  interest  is  indicated  in  figure  2.  The  m-th  local 
reference,  rm(t),  can  be  taken  equal  to  s  (t);  a  clipped  version, 
such  as  sgn{sm(t)};  or  some  other  distorted  version,  such  as 
g{sm(t)},  if  desired. 


K  inputs  Clip  Sum  on  Multiply  Sum  on  Output 

Space  References  Time  Variables 


Figure  2.  DIMUS  Processor 


T 

The  output  variables  Z  s  [z^,...,zM]  in  figure  2  are  subject 
to  additional  processing  before  a  decision  is  made  as  to  which 
particular  signal  was  transmitted.  Just  as  for  the  earlier 
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linear  processor  (3),  output  vector  Z  is  subject  to  a  linear 
transformation  into  decision  vector  Z  =  A  Z  +  B,  where  A  and  B 
are  known  matrices.  Identity  decisions  are  made  upon  Z.  This  is 
not  necessarily  the  optimum  processing  to  conduct  upon  DIMUS 
output  Z;  rather,  guided  by  the  linear  processor,  it  is  adopted 
as  a  reasonable  procedure  to  employ  for  the  DIMUS  processor  as 
well . 

When  the  product  of  the  number  of  channels  and  the  number  of 
time  samples  in  figure  2,  namely,  KT,  is  large,  each  output 
variable  z  is  the  sum  of  a  large  number  of  random  variables. 

In  this  case,  the  random  output  vector  Z  tends  toward  a  set  of 
joint  Gaussian  random  variables;  for  a  justification,  see  the 
extensive  simulation  using  K  *  256  and  T  =  4096  in  reference  1, 
page  20. 


The  above  observation  allows  the  analytic  investigation  of 
the  DIMUS  processor  to  concentrate  on  finding  the  mean  vector  and 
the  covariance  matrix  of  output  Z.  Then,  the  corresponding  mean 
and  covariance  statistics  for  decision  vector  Z  are  easily 
determined  by  means  of  equation  (4),  at  which  point  a  low-order 
M-ary  simulation  of  classification  performance  can  be  conducted. 
The  comments  made  earlier  in  the  Simulation  subsection  are 
directly  relevant  here. 
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ANALYSIS 


The  m-th  output  in  figure  2  is  given  by 


T  K 

zm  "  ZZH  rm(t)  YZ  yk(t)  for  1  <  M  ,  (10) 

m  t-1  m  k=l  K 


where 


yk(t)  =  sgn{sm(t)  +  nk(t)}  for  1  <  k  <  K  .  (11) 


Using  the  zero-mean  Gaussian  character  of  input  noises  (nk(t)} 
shows  the  mean  of  z  to  be  given  by 


T  K 

E{z  }  =  YU  rm(t)  YU  E(y,  (t)}  for  1  <  m  <  m  ,  (12) 

m  t-1  k=l  K 


with 


E{yk<t)}  =  E{sgn{sm(t)  +  nk(t)}}  = 

2 

=  J  dx  sgn{sro(t)  +  x}  ^2ncrkj-!?  exp(-  = 

—  2  a, 


-  2  ^(sm(t)/crkj  -  1  for  1  <  k  <  K  , 


(13) 


2  2 

where  ck  is  the  noise  variance  E{nk(t)}  =  Rkk(0)  in  the  k-th 
channel,  and  normal  probability  function  $  and  PDF  <j>  are  given  by 


*( 


x)  s  J  dt  (2n)  ^  exp(-t^/2)  =  J  dt  <j>(t)  . 


(14) 


Substitution  of  equation  (13)  into  equation  (12)  yields  mean 
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E‘%) 


T 

rz 

t-i 


K 

c 

k-l 


|^2  $  {&m(  t  )/crkj -1  ]  for  1  <  m  <  M  .  (15) 


This  result  requires  a  double  sura  over  KT  terms,  just  as  the 
output  zm  itself  does  in  equation  (10).  In  the  special  case 
where  all  the  input  noise  standard  deviations  crk  are  equal  to  a, 
equation  (15)  simplifies  to 


E  { z  }  = 
'  m 


T 

K  C  rm(t)  [2*(sm(t)/<r)-l]  for  1  <  m  < 


M 


(16) 


Comparison  of  equation  (16)  with  the  analogous  result  (6)  for 
the  linear  processor  reveals  that  sm(t)  is  now  replaced  by  a 
memoryless,  nonlinearly  distorted  version.  The  function  2$(x)-l 

\s 

is  linear  for  small  x,  namely,  (2/n)  2  x;  however,  it  saturates  at 
±1  as  x  ->  ±<°. 


In  order  to  determine  the  covariance  matrix  of  DIMUS  output 


vector  Z,  consider  the  crosscorrelation  of  output  variables  zm 
and  z*  from  equation  (10): 


Efzm  ZX>  -  c  r  (t)  rx(u)  C  E(y  (t)  y^u))  .  (17) 

t , u=l  k , ] =1  J 


From  equation  (11),  the  required  ensemble  average  follows  as 


E { yk ( t )  y j ( u ) }  =  E{sgn{sm(t)  +  nk(t)}  sgn{sm(u)  +  n..(u)}}  = 

=  JJ  dx  dy  sgn{sm(t)  +  x}  sgn{sm(u)  +  y}  pk_.  (x,y,t-u)  ,  (18) 


where  the  joint  second-order  stationary  Gaussian  noise  PDF  pj.  is 
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Covariance  coefficient  pkj(x)  =  j  (  t)/(  <rk  a  j  )  .  Evaluation  of 
double  integral  (18)  is  conducted  in  appendix  B,  with  infinite 
expansion  (B-7)  as  the  end  result. 


The  covariance  of  interest  is  given  by  equation  (B-8)  as 
Cov{yk(  t)  ,y_.(u) }  =  c(pkj(t-u)  »sm(  t  )/ffk,sm(u)/(j;.]  (20) 


for  1  <  k,j  <  K,  where  function 


n+1 


C(p,x,y)  .  |  exp(-  X  -y-*  )  C  ( n+! )  1  HenU)  Henty)  •  (21) 


A  BASIC  program  for  this  form  of  function  C(p,x,y)  is  provided  in 
appendix  B.  Additional  forms  for  covariance  Cov{yk( t ) ,y^  (  u) } ,  as 
well  as  for  the  limiting  cases  of  p  =  ±1,  are  also  presented  in 
this  appendix. 


The  covariance  of  DIMUS  outputs  zm  and  z^  now  follows  from 
equations  (17)  and  (20)  as 


Cov‘Vz!l)  ’  C  rm<t>  EZ  c°v{yk(t),y.(u))  - 

t  ,u=l  k , 3=1  J 

T  K 

=  C  r  (U  rx(u)  C  c(pkj(t-U),sIn(t)/Pk,snl(u)/aj)  (22) 

l  ^  ll  —  i  K  f  ]  “J. 


for  1  <  m,X  <  M. 
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Equations  (15)  and  (22)  furnish  the  mean  and  covariance 

information  required  to  statistically  characterize  the  DIMUS 

output  random  variables  in  figure  2,  at  least  for  large  KT. 

However,  when  that  is  true,  equation  (22),  which  requires  that 
2 

(KT)  terms  be  calculated  and  summed  for  each  pair  m,X,  will 
entail  a  considerable  amount  of  calculations. 


If  the  expression  for  function  C  in  equation  (21)  is 
substituted  into  equation  (22),  and  if  all  the  noise  standard 
deviations  are  equal  to  a,  the  covariance  can  be  expressed  as 


OO  "P 

c°v(zm,zx}  =  (T^rrr  , rm< t >  Vu>  exp(- 

n=u 


a2 ( t )  +  b2 ( u ) ' 


x  Hen(a(t))  Hen(b(u))  Fn(t-u)  for  1  <  m,X  <  M  ,  (23) 


where 

a  ( t )  ■  sm(t)/c  ,  b(u)  s  sm(u)/<r 


F„(t) 


K 

-  L2 

k,  j-1 


n+1 i 
Pkj  (T) 


(24) 


If  the  double  sum  on  k , j  in  equation  (24)  can  be  evaluated  in 
closed  form  (for  some  noise  spectra),  expansion  (23)  may  afford  a 
quicker  method  of  evaluation  than  equation  (22).  Recursions  on 
the  Hermite  functions  in  expansion  (23)  could  still  be  used  to 
advantage;  however,  most  likely,  there  will  not  be  a  recursion 
for  sequence  {Fn(x)}. 


14 


EXECUTION  TIME 


Some  sample  execution  times  for  the  evaluation  of  the 
covariance  in  equation  (22)  using  a  Hewlett-Packard  9000  computer 
are  listed  in  table  1.  These  times  are  for  a  single  value  of  the 
m,X  pair,  and  it  is  presumed  that  the  rm(t),  r^(u),  p^^Ct-u),  and 
s  (t)/cr,  arrays  have  already  been  computed  and  filled  in.  On  the 
other  hand,  no  use  of  any  inherent  symmetries  have  been  taken 
advantage  of,  which  would  have  reduced  the  execution  time; 
rather,  equation  (22)  was  evaluated  directly  in  the  form  given. 


Table  1.  Sample  Execution  Times 


K 

T 

Time  (sec) 

10 

12 

41.1 

10 

24 

166.7 

20 

24 

656.9 

For  a  CRAY  computer,  these  times  would  be  reduced  by  almost  a 

factor  of  1000,  requiring  about  1  second  for  K  =  20,  T  =  24.  It 

can  be  observed  that  the  execution  times  behave  proportionally  to 
2  2 

K  and  T  .  Thus,  K  and/or  T  could  be  significantly  increased  for 
a  comparative  study  of  linear  and  DIMUS  beamformers.  However, 
the  earlier  case  of  K  =  256,  T  =  4096  (reference  1)  appears  to 
require  excessive  computation  time,  even  for  a  CRAY  computer. 


The  execution  time  is  heavily  dependent  on  the  sizes  of  the 
covariance  coefficients  {  p j  ( t— u ) }  encountered  in  equation  (22), 
as  well  as  on  the  error  tolerance  (Tol)  used  in  the  program  in 
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appendix  B.  The  timing  results  above  are  for  a  uniform 
distribution  of  p  values  over  +0.8  and  Tol  =  l.E-8.  If 
p  values  are  uniformly  distributed  over  +0.5  and  Tol  is 
to  l.E-6,  the  execution  time  is  reduced  by  half. 


the 

increased 
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APPENDIX  A  -  GENERATION  OF  CORRELATED  RANDOM  VARIABLES 
WITH  SPECIFIED  MEAN  AND  COVARIANCE  MATRICES 

T 

Let  U  =  [u^,...,uM]  be  a  set  of  M  zero-mean,  unit-variance 
uncorrelated  random  variables.  Form  the  linear  transformation 

V  =  A  U  +  B  ,  ( A— 1 ) 

where  matrix  A  is  MxM  and  matrix  B  is  Mxl .  Then,  random  Mxl 
vector  V  has  mean  and  covariance  matrices  given  by 

E{ V}  =  V  =  B  , 

Cov{  V}  =  A  U  UT  AT  =  A  AT  ,  ( A— 2  ) 

respectively,  where  an  overbar  denotes  an  ensemble  average. 

If  the  mean  and  covariance  matrices  of  vector  V  are  specified 

as  Uv  and  Cv,  respectively,  the  two  matrices  in  equation  (A-l) 

T 

must  be  taken  according  to  B  =  Uv  and  A  A  =  Cy.  That  is,  A  is  a 
square  root  matrix  of  Cv . 

If,  in  addition,  all  the  random  variables  (u^)  in  equation 
(A-l)  are  Gaussian,  then  random  variables  V  are  joint  Gaussian, 
because  equation  (A-l)  is  a  linear  transformation.  In  this  case, 
equation  (A-2)  yields  a  complete  statistical  description  of 
random  variables  V. 


A-l/( A-2  blank) 


APPENDIX  B  -  EVALUATION  OF  DOUBLE  INTEGRAL  (18) 


The  desired  crosscorrelation  is  obtained  by  substituting 
joint  PDF  (19)  into  equation  (18).  Letting 

a  =  s  ( t )  /  o,  ,  b  =  s  ( u )  /  c .  ,  x '  =  x/u,  ,  y'  =  y/a.  ,  (B-l) 


and  then  replacing  x'  and  y'  by  x  and  y,  respectively,  for 
notational  simplicity  results  in 

E{yk( t)  Yj (u) }  =  JJ  dx  dy  sgn{a  +  x}  sgn{b  +  y}  x 


2.2 


2nd  -  p  )' 


L  .  p  ,  p  (t-u)  .  (B-2) 

V  2(1  -  pZ)  J  K3 


At  this  point,  the  following  expansion  (reference  2,  chapter 
26)  is  employed,  which  separates  the  dependence  on  x,y,p: 


2.  2 


2.2.  n 


„p(_  *  +*  -  exp  (-  5_p)  E  He  (x )  Hen(y) 

(1  -  t>  K  2(1  -  o  )  1  '  ’  n-0 


Substitution  in  double  integral  (B-2)  yields 


( B— 3 ) 


Uyk<t)  yj(u>)  -  C  t  En(a)  En(b) 


(B-4) 


where 


En(a)  s  J  dx  sgn{a  +  x}  (2n)  ^  exp(-x2/2)  Hen(x)  .  (B-5) 

There  follows,  upon  use  of  equation  (14),  Eg(a)  =  2$(a)-l.  By 
use  of  reference  2,  equation  22.11.8,  the  remaining  functions  are 


En ( a )  =  exp(-a2/2)  (-l)n_1  Hen_1(a)  for  n  >  1  .  (B-6) 


B-l 


Collecting  these  results,  equation  (B-4)  becomes 


E{yk(t)  Yj (u) }  =  [ 2$(a)-l ]  [2$(b)-l]  + 

+  i  exP(-  2  b  ')  C  Hen-l(a)  Hen-l(b)  •  (B-7) 

n=i 

(This  expansion  is  essentially  the  same  as  reference  2,  equation 
26.3.29.)  When  the  identity  of  a  and  b  in  equation  (B-l)  is 
recalled,  the  leading  term  in  equation  (B-7)  is  recognized,  by 
reference  to  equation  (13),  as  the  product  of  the  mean  values  of 
y,  (t)  and  y.(u).  Then,  shifting  integer  n  by  1,  the  covariance 

K  j 

follows  from  equation  (B-7)  as 

Cov{yk(t)  ,y_.  (u) }  ■  c(pkj(t-u),sm(t)/<rk,sm(u)/«Xj]  (B-8) 

for  1  <  k,j  <  K,  where  function 

°°  n+1 

C  Tferrr  Hen(x)  Hen(7)  •  (B-9> 

The  Hermite  polynomial  Hen(x)  grows  approximately  as  (n!)2 
for  increasing  n;  see  reference  3,  equation  8.22.8.  Therefore, 
direct  recursions  on  each  of  the  three  multiplicative  terms  in 
the  sum  of  equation  (B-9)  will  likely  encounter  underflow  and 
overflow  before  sufficient  convergence  is  achieved.  In  order  to 
alleviate  these  effects,  a  factor  T(n/2+l)  2n ^  was  inserted  in 
the  denominators  of  both  Hermite  functions,  and  the  square  of 
that  factor  was  applied  to  (multiplied  by)  the  first  term.  This 
factor  increases  in  a  similar  fashion  to  ( n l ) 2  for  large  n,  but 
can  be  evaluated  by  a  recursive  procedure  without  square  roots. 


C( p,x ,y)  s  ~ 


M-  H1-) 


B— 2 


With  the  definitions 


-  «p<-uV)/«>  .  *n  -  - 


b  r( ( n+1 )/2 )  z  Hen(z)  E 

n  V2  r( n/2+1 )  '  n 


r ( n/2+1 )  2 


n/2 


,  z  =  x  or  y  ,  ( B— 1 0 ) 


expansion  ( B— 9 )  can  be  expressed  as 


C  (  p ,  x  ,  y )  -  TZ  An  R*  B?n  , 
n=0 


( B— 11 ) 


with  starting  values  and  recursions  (for  n  >  2) 


A0  -  t’  ■  A1  -  Ip2  '  A-  -  A 


2  n 


n  n-2  K  n+1  ' 


-fir  R  -  f— 1 2  B  -  B  . 

UJ  '  B1  “  UJ  '  Bn  "  Bn-2  n  ' 


R0  -  E  ■  R1  -  E  z  (i f  ■  Rn  Rn-1  z  Bn  '  Rn-2  ^  •  (B-12> 


The  last  recursion  was  derived  from  reference  2,  equation 
22.7.14,  and  definition  ( B— 1 0 )  above. 


An  efficient  BASIC  program  for  C(p,x,y)  that  incorporates 
these  modified  functions  is  listed  in  figure  B-l;  it  employs  only 
the  one  exponential  E  in  equation  (B-10).  The  tolerance  (Tol) 
and  maximum  number  of  terms  utilized  (Num)  are  under  the  user's 
control.  The  tolerance  requirement  is  imposed  on  the  sum  of  the 
magnitudes  of  the  latest  two  terms  in  the  series,  relative  to  the 
current  sum  at  that  point. 


B-3 


10  DEF  FNC ( P , X , Y ) 

20  T=X*X+Y*Y 

30  IF  ( T>2  0  0  . )  THEN  RETURN  0. 

40  Tol=l . E-8  !  SPECIFIED  RELATIVE  ERROR 

50  Num=100  0  !  MAXIMUM  NUMBER  OF  TERMS 

60  E=EXP ( — . 25*T ) 

70  P2=P*P 

80  A2=. 63661977236758134*P 

90  Al= . 5*P2 

100  B2=l .2533141373155003 

110  Bl=l./B2 

120  R2x=E 

130  Rlx=E*X*Bl 

140  R2y=E 

150  Rly=E* Y*Bl 

160  Fl= . 5 

170  Tl  =  l . El  0  0 

180  S=A2*R2x*R2y * (l.+.5*P*X*Y) 

190  FOR  N=2  TO  Num 
200  F=N/(N+1) 

210  A=A2*P2*F 

220  B=B2*Fl 

230  Rx=Rlx*X*B— R2x*Fl 

240  Ry=Rly * Y*B-R2y*Fl 

250  T=A*Rx*Ry 

260  S=S+T 

270  IF  ( Tl+ABS ( T )  ) <= ( Tol *AB S ( S )  )  THEN  410 

280  A2=Al 

290  Al=A 

300  B2=Bl 

310  Bl=B 

320  R2x=Rlx 

330  Rlx=Rx 

340  R2y=Rly 

350  Rly=Ry 

360  Fl=F 

370  T1=ABS(T) 

380  NEXT  N 

390  PRINT  "TOLERANCE  NOT  SATISFIED" 

400  PAUSE 
410  RETURN  S 
420  FNEND 


Figure  B-l.  Program  for  Function  C(p,x,y)  in  Equation  (B— 9) 


B-4 


With  the  aid  of  equations  (13),  (14),  (B-l),  and  (B-2),  the 
covariance  can  be  manipulated  into  the  single  integral  form 

00 

Cov{yk(t),yj(u)}  =  4  J  dx  +(x)  [*(b)  -  *(  b  ~  p*  J]  ;  (B-13) 

3  ^ 

the  alternative  form  with  a  and  b  switched  is  also  allowed. 

For  p  =  ±1,  none  of  the  above  forms  are  directly  usable. 
Instead,  a  direct  evaluation  of  the  covariance  in  these  two 
special  cases  yields 

'4  $(a)  $(-b)  for  a  <  b' 

Cov  =  <  .  for  p  =  +1  ( B-l 4 ) 

.4  $(-a)  ${b)  for  a  >  b, 

and 

-4  $(a)  $(b)  for  a+b  <  O' 

Cov  =  «  .  for  p  =  -1  .  ( B-15 ) 

-4  $(-a)  $(-b)  for  a+b  >  0, 


B— 5/( B-6  blank) 
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